Spectra of Random Contractions and Scattering Theory for 
Discrete-Time Systems 



Yan V. Fyodorov^f, § and Hans-Jiirgen Sommers§ 
Department of Mathematical Sciences, Brunei University, Cleveland Road, Uxbridge, UB8 3PH, UK 
%Fachbereich Phystk, Universitat-GH Essen, D-45117 Essen, Germany 
(February 8, 2008) 

O ' 

Random contractions (sub-unitary random matrices) appear naturally when con- 
sidering quantized chaotic maps within a general theory of open linear stationary 
systems with discrete time. We analyze statistical properties of complex eigenval- 
ues of generic N x N random matrices A of such a type, corresponding to systems 
with broken time-reversal invariance. Deviations from unitarity are characterized 
by rank M < N and a set of eigenvalues < Ti < l,i = 1, M of the matrix 
on ; f = i - i 1 " A. We solve the problem completely by deriving the joint probability 

density of ./V complex eigenvalues and calculating all n— point correlation functions. 
In the limit N » M,n the correlation functions acquire the universal form found 
earlier for weakly non-Hermitian random matrices. 
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The theory of wave scattering can be looked at as an integral part of the general theory of 
linear dynamic open systems in terms of the input-output approach. These ideas and relations 
were developed in system theory and engineering mathematics many years ago, see papers [0-j3| 
I ' and references therein. Unfortunately, that development went almost unnoticed by the majority 
' of physicists working in the theory of chaotic quantum scattering and related phenomena, see [^J5| 
and references therein. For this reason we feel it could be useful to recall some basic facts of the 

IT) , 

q°7** 3 input-output approach in such a context. 

An Open Linear System is characterized by three Hilbert spaces: the space Eq of internal states 
$ G Eq and two spaces E± of incoming (-) and outgoing (+) signals or waves also called input 
and output spaces, made of vectors <p± G E±. Acting in these three spaces are four operators, or 
matrices: a) the so-called fundamental operator A which maps any vector from internal space Eq 
onto some vector from the same space Eq b) two operators Wi,%, with W\ mapping incoming states 
onto an internal state and W% mapping internal states onto outgoing states and c) an operator Sq 
acting from E- to E + . 

We will be interested in describing the dynamics ^f(t) of an internal state with time t provided 
we know the state at initial instant t = and the system is subject to a given input signal (j)-(t). 
In what follows we consider only the case of the so-called stationary (or time-invariant) systems 
when the operators are assumed to be time- independent. Let us begin with the case of continuous- 
time description. The requirements of linearity, causality and stationarity lead to a system of two 
dynamical equations: 

i£* = i*(t)+Wi^_(t) 

+ (t)=£ O 0-(t)+*W2*(t) 

Interpretation of these equations depends on the nature of the state vector \I> as well as of the 
vectors <f>± and is different in different applications. In the context of quantum mechanics one 
relates the scalar product with the probability to find a particle inside the "inner" region at 
time t, whereas 4>±<l>± stays for probability currents flowing in and out of the region of internal 
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states (the number of particles coming or leaving the inner domain per unit time). The condition 
of particle conservation then reads as: 

AvE-t^ = _ (2) 

It is easy to verify that Eq. (|2j) is compatible with the dynamics Eq.(|l|) only provided the operators 
satisfy the following relations: 

A ] - A = iWW^ , S f S a = i andW^ = W 2 = -S W\ 

which shows, in particular, that A can be written as A = H — ~WW\ with a Hermitian H = H<. 

The meaning of H is transparent: it governs the evolution i-j^ty = H^(t) of an inner state \& 
when the coupling W between the inner space and input /output spaces is absent. As such, it is 
just the Hamiltonian describing the closed inner region. The fundamental operator A then has 
a natural interpretation of the effective non-self adjoint Hamiltionian describing the decay of the 
probability from the inner region at zero input signal: <fi-(t) = for any t > 0. If, however, the 
input signal is given in the Fourier-domain by 0_(u>), the output signal is related to it by: 

(f>+{uj) = \§(w)So] <M«) , S{uj) = i - iW* — - W (3) 
L J col — A 

where we assumed ^S(t = 0) = 0. The unitary matrix S(u>) is known in the mathematical literature 
as the characteristic matrix-function of the non-Hermitian operator A. In the present context it is 
just the scattering matrix whose unitarity is guaranteed by the conservation law Eq(||). 

The contact with the theory of chaotic scattering is now apparent: the expression Eq.(^|) was 
frequently used in the physical literature [^J^] as a starting point for extracting universal properties 
of the scattering matrix for a quantum chaotic system within the so-called random matrix approach. 
The main idea underlying such an approach is to replace the actual Hamiltonian H by & large 
random matrix and to calculate the ensuing statistics of the scattering matrix. The physical 
arguments in favor of such a replacement can be found in the cited literature. 

In particular, most recently the statistical properties of complex eigenvalues of the operator A 
as well as related quantities were studied in much detail [|| |sj . Those eigenvalues are poles of the 
scattering matrix and have the physical interpretation of resonances - long-lived intermediate states 
to which discrete energy levels of the closed system are transformed due to coupling to continua. 

In the theory presented above the time t was a continuous parameter. On the other hand, a 
very useful instrument in the analysis of classical Hamiltonian systems with chaotic dynamics are 
the so-called area-preserving chaotic maps . They appear naturally either as a mapping of the 
Poincare section onto itself, or as result of a stroboscopic description of Hamiltonians which are 
periodic functions of time. Their quantum mechanical analogues are unitary operators which act 
on Hilbert spaces of finite large dimension N . They are often referred to as evolution, scattering 
or Floquet operators, depending on the physical context where they are used. Their eigenvalues 
consist of N points on the unit circle (eigenphases). Numerical studies of various classically chaotic 
systems suggest that the eigenphases conform statistically quite accurately the results obtained for 
unitary random matrices of a particular symmetry (Dyson circular ensembles). 

Let us now imagine that a system represented by a chaotic map (" inner world" ) is embedded in 
a larger physical system (" outer world" ) in such a way that it describes particles which can come 
inside the region of chaotic motion and leave it after some time. Models of such type appeared, for 
example, in JTo[ ] where a kind of scattering theory for " open quantum maps" was developed based 
on a variant of Lipmann-Schwinger equation. 
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On the other hand, in the general system theory |2|3j] dynamical systems with discrete time 
are considered as frequently as those with continuous time. For linear systems a "stroboscopic" 
dynamics is just a linear map (0_(n); ^>(n)) — > (</> + (n); ^(n + 1)) which can be generally written 



as: 






= V > £ 7 (4) 



Again, we would like to consider a conservative system, and the discrete-time analogue of Eq. (g) 
is: 

*t( n + l)*(n + 1) - * f (n)*(n) = <p f _(n)(f>-(n) - ^ + {n)(j) + (n) 

which amounts to unitarity of the matrix V in Eq.(0). In view of such a unitarity V of the type 
entering Eq.(^) can always be parametrized as (cf. PM): 



\o v x J y ft vi-fTf y ^ o v 2 J w 

where the matrices and 11^2 are unitary and f is a rectangular N x M diagonal matrix with 
the entries = Sij Tj , 1 < i < N , 1 < j < M < Tj < 1 . 

Actually, it is frequently convenient to redefine input, output and internal state as: 4>—{n) — > 
<j>—(n), 0+ — * v\4>+{n) and \&(n) — > U2 , 5(n) which amounts just to choosing an appropriate 
basis in the corresponding spaces. The transformations bring the matrix V to a somewhat simplier 
form: 



where u = u\u\. Such a form suggests a clear interpretation of the constituents of the model. 
Indeed, for f = the dynamics of the system amounts to: \&(n + 1) = u^(n). We therefore 
identify u as a unitary evolution operator of the " closed" inner state domain decoupled both from 
input and output spaces. Correspondingly, f ^ just provides a coupling that makes the system 
open and converts the fundamental operator A — uV 1 — f f t to a contraction: 1 — A' A = tt^ > 0. 
As a result, the equation fy(n + 1) = Aty(n) describes an irreversible decay of any initial state 
^(0) for zero input (f>—(n) = 0, whereas for a nonzero input and ^(0) = the Fourier-transforms 
4>±(oj) — X^^Lo e maJ 0±(n) are related by a unitary scattering matrix S(u>) given by: 



S(u) = >/l-ftf - f f r &f (7) 

Assuming further that the motion outside the inner region is regular, we should be able to 
describe generic features of open quantized chaotic maps choosing the matrix u to be a member 

■ 2 

, with the bar standing for the 



of a Dyson circular ensemble. Then one finds: f'f = 1 — S(uj) 

averaging of §(w) in Eq.(@) over u. Comparing this result with MM we see that the M eigenvalues 
< Ti < 1 of the M x M matrix ftf play the role of the so-called transmission coefficients and 
describe a particular way the chaotic region is coupled to the outer world. 

In fact, this line of reasoning is motivated by recent papers Jl^,[l3| • The authors of |l2| considered 
the Floquet description of a Bloch particle in a constant force and periodic driving. After some 
approximations the evolution of the system is described by a mapping: c„ + i = Fc„, where the 
unitary Floquet operator F = SU is the product of a unitary "M-shift" S : Ski — Si.k-M-, l,k = 
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— oo, oo and a unitary matrix U . The latter is effectively of the form U = diag(di, u, (I2), where 
^1,2 are (semi)infinite diagonal matrices and u can be taken from the ensemble of random N x N 
unitary matrices. 

One can check that such a dynamics can be easily brought to the standard Eqs.(0,H) with the 



fundamental operator being an N x N random matrix of the form A = yl — ftfu, and all M 
diagonal elements of the N x M matrix t are equal to unity. Actually, the original paper |l2| 
employed a slightly different but equivalent construction dealing with an "enlarged" internal space 
of the dimension N + M. We prefer to follow the general scheme because of its conceptual clarity. 

Direct inspection immediately shows that the non-vanishing eigenvalues of the fundamental 
operator A as above coincide with those of a (N — M) x (N — M) subblock of the random unitary 
matrix u. Complex eigenvalues of such "truncations" of random unitary matrices were studied in 
much detail by the authors of a recent insightful paper fll3|] . They managed to study eigenvalue 
correlations analytically for arbitrary N, M. In particular, they found that in the limit N — > 00 
for fixed M these correlation functions practically coincide |l2| with those obtained earlier ||^| for 
operators of the form A = H — jWW^ occuring in the theory of open systems with continuous-time 
dynamics. 

Such a remarkable universality, though not completely unexpected, deserves to be studied in more 
detail. In fact, truncated unitary matrices represent only a particular case of random contractions 
A. Actually, some statistical properties of general subunitary matrices were under investigation 
recently as a model of scattering matrix for systems with absorption, see [[uj. 

The particular case of rank-one deviations from unitarity is the simplest one to investigate and 



was considered in a recent preprint |15|. However, generalization to arbitrary M along the lines of 
flEf seemed to be problematic. The main goal of the present paper is to suggest a regular way of 
studying the specta of random contractions for a given deviation from unitarity. 

The ensemble of general N x N random contractions A = iiy/l — ff t describing the chaotic map 
with broken time-reversal symmetry can be described by the following probability measure in the 
matrix space: 

V{A)dA cx 6(A*A - G)dA , G = l- f (8) 

where dA = rj f . dAijdA*^ and we assumed that the unitary matrix u is taken from the Dyson 
circular unitary ensemble. The NxN matrix f f t = 1 — G > is natural to call the deviation matrix 
and we denote it T. It has M nonzero eigenvalues coinciding with the transmission coefficients T a 
introduced above. The particular choice Ti<M — 1, 2i>M = corresponds to the case considered 
in [T^ ] . In what follows we assume all Tj < 1 , but the resulting expressions turn out to be valid in 
the limiting case T, = 1 as well. 

Our first step is, following |^ Jl^ , |l^| , introduce the Schur decomposition A = U(Z + R)U^ of 
the matrix A in terms of a unitary U, diagonal matrix of the eigenvalues Z and a lower tri- 
angular R. One can satisfy oneself, that the eigenvalues z%,...,Zn are generically not degen- 
erate, provided all Tj < 1. Then, the measure written in terms of new variables is given by 
dA = I A({z})\ 2 dRdZdfi(U), where the first factor is just the Vandcrmonde determinant of eigen- 
values Zi and dfj,(U) is the invariant measure on the unitary group. The joint probability density 
of complex eigenvalues is then given by: 

V({z}) cx |A({z})| 2 [ d^{U)dR 8 Uz + R){Z + i?)t - U^GU^j (9) 



The integration over R can be performed with some manipulations using its triangularity (some 
useful hints can be found in |13|). As the result, we arrive at: 
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r N 

P({z})<x |A({z})| 2 / d f i(U)Y[s(\z 1 \ 2 ...\z l \ 2 -det [l-TUP^]) (10) 

where P = diag (1, ...1,0, 0), with first I entries being equal to unity, is a projector. 

To perform the remaining integration over the unitary group was mentioned as the main technical 
problem in jy| , and proved to be more difficult than the corresponding procedure for non-Hermitian 
matrices, see Q. Here we outline the main steps and present some intermediate expressions 
relegating details of the calculation to a more extended publication. First, one considers the 
columns of N x N unitary matrices as N— component of (mutually orthogonal) vectors aj, I = 
1, iV and introduces N x N matrices of rank I: Qi = a « ® a l i so ^ na t UPiU^ = Qi. Writing 

down the corresponding constraints in a form of 8— functions, one can represent the expression 
Eq.(||) in the form: 

r Z^ -1 \ N N r 

V({z})<x |A({z})| 2 / J] dQi Jp (M 2 -M 2 - det (! - TQij J[ I dadatS (4 - Qi-i - a ® at) (11) 

\ 1=1 / 1=1 i=l 

where the matrices Qi are considered to be unconstrained N x N Hermitian. We also used the 
orthonormality condition Qn = 1 as well as the convention Qq = 0. 

Due to the fact that only M out of N eigenvalues of the matrix T are non-zero, both the matrices 
Qi and the vectors a can be effectively taken to be of the size M ( it amounts to changing the 



unspecified normalisation constant in Eq.(ll)) and redefine the matrix T as T — diag(Ti, Tm) = 
ftf. Then it is convenient to change: T^-^QiT 1 / 2 — > Qi and separate integration over eigenvalues 
and eigenvectors of matrices Qi. The latter can be performed in a recursive way I — > I + 1, with 
the multiple use of the Itzykson-Zuber-Harish-Chandra [|l6| formula. After quite an elaborate 
manipulation, one finally arrives at the following representation: 



where we defined the diagonal matrices of size M as: t — diag(ri, ...,tm) , A = diag(Ai, Am) 
and used the notations: r c = In (1 — T c ) and 

/ (a,A)=z^ 1 ^ 6 . (13) 

The distribution Eq.(|l2|) is written for |zfe| 2 < 1 for any k — 1,...,N and vanishes otherwise. 
The remarkable feature of such a distribution is that it allows for calculation of all n— point corre- 
lation functions for arbitary N, n, M with help of the method of orthogonal polynomials. Again, 
the particular case M = 1 (lj| is quite instructive and can be recommended to follow first for 
understanding of the general formulae outlined below. 

To this end, we write 



N N 

iA({z»i 2 n m 2 > *) = n det 

k=l k=l 



N 



h N n (X) 



(14) 

i,j = l,...N 



where the constants N n (A) are provided by the orthonormality condition: 

d 2 ^ m - 1 (^)«- 1 /(ln|^| 2 ,A)=<5 m)n iV n (A) , (15) 



p<i 



G 



which yields after a simple calculation N n (X) = nY[ c= i ( n +i\ ) • 

Now, by applying the standard machinery of orthogonal polynomials fll7| one can find the cor 
relation function: 

Nl 



Rn (-2*1 j **■) %7i) 



[N — n)\ 



d z n+1 ...d z N T>{{z}) 



N 



as given by: 

R n (zi, z n )<xV f d\e~ lTTf ~ x Y[ N k (X) det \k(z u Zj ; A) 

^ k=l 

where the kernel K is defined as: 

1 N 

K( Zl ,z 2] A) = - V det (iX + n){ Zl zZ) n - l f{\n |z 2 | 2 , A) 

TT ' * 



(16) 



(17) 



(18) 



and the differential operator T> is just the expression in front of the A— integral in Eq.([L2|). 

In principle, all A— integrations in the equation Eq. (|l7j) can be performed explicitly and the 
resulting formulae provide the desired general solution of the problem. However, for arbitary 
N, M, n the results obtained in that way are still quite cumbersome. We present below as an 
example the lowest correlation function R\{z), which is just the mean eigenvalue density inside 
the unit circle \z\ < 1. It can be calculated from the following recursive relation connecting the 
density for M and M — 1 open channels: 



r[ M) (z) = RY vl ~ 1 >(z) + 



(M-l). 



1 d 



(M-l), 



where 



t[ m) {t c m 2 }=y: 



i - 



M 2 



N-l 



1+T C ) 



4 M -V{T c ;\z\ 2 } 



I |2JV roc M-l , 1 

(N- 1)! J l\\T c 



1 d 



t d\i 



1- z 



i2N 



1- Z 



(19) 



(20) 



(21) 



For the case of all equivalent channels, i.e. when all the transmission coefficients T c are equal: 
T c = T, such a recursive relation can be represented in a more compact form: 



R [ M \z) = r[ M ~ 1 \z) + 



1 



d 



(M-l)! \dt 
where the generating function lZ{t, z) is given by: 



M-l 



K{t,z)\ t= o , 



K(t,z) 



ir d\z\ 2 £ — rj 



and 



C = l + (t - 1 



T 



(22) 



(23) 



(24) 



All the equations above are valid for arbitrary N > M, n. In the theory of quantum chaotic 
scattering we, however, expect a kind of universality in the semiclassical limit. Translated to the 
random matrix language such a limit corresponds to N — > oo at fixed n, M . Still, extracting the 
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asymptotic behaviour of the correlation function R n (z\ 1 ...,z n ) from Eq.([L7|) in such a limit is not 
a completely straightforward task. A useful trick is to notice that Eq.(n7n can be rewritten as: 



M 



Rk(zi, ...,z n ) oc 

q 1= l,...,q n = l 

N 



F qu ..., qn ({Tc,z}) = B{T C ] det 



On 



9tm 



+ k ) (z l z*) k - 1 



(25) 
(26) 



x 



n(=- 



ci<c 2 

where we introduced the notation: 



dr Cl 8t C2 



U[dX c 



exp —i\ c T c 



l,...,n 
n TTm 



a/ 



£{T C } 



rpM—N 

n ci<C2 (T C1 -T C2 )1 = 1 (1-T C ) 



Introducing now the auxiliary differential operator = IIj=i Flf/g (aT - — af - ) an< ^ con ~ 

sidering its action upon the ratio F qii , ^ qri /B{T C } one can satisfy oneself that in the limit N 3> M, n 



the leading contribution to Fq lt ...,q n is given by: 



M 



(1 - T c ) / f c 



N-M 



M 



c=l 



where the kernel is given by 

K(z l ,z J ;{f c })=J2U 



nn 

j=l s^ qj 



1 1 



det 



K(zi,Zj-,{f c }) 



N M 



k=l c=l 
M 



(N — M) - . Tfc + /■■ - 1 



n 



(N-M) 



1-Te 

T, 



1-sb i: 



* \ /c — 1 



(28) 



where we used the notation: T c = 1 — exp (t c — Y^j=i ^qj.c m \ z j\ 2 

Further simplifications occur after taking into account that eigenvalues Zi are, in fact, concen- 
trated typically at distances of order of 1/N from the unit circle. Then it is natural to introduce 
new variables yi,(f>i according to = (1 — yi/N)e b ^ i and consider yi to be of the order of unity 
when N — > oo. First of all, one immediately finds that: 



M 



lim TT ( ^ 



N-M 



exp 



c=l 



-»£-- 1 - r - 



3 = 1 



T„ 



(29) 



As to the phases fa, we expect their typical separation scaling as: fa — fa — 0(1 /N). Now it is 
straightforward to perform explicitly the limit N — > oo in Eq.(|8|). Combining all factors together, 
one brings the correlation function Eq. (p5|) to the final form: 



e ~9qVk 



71 M 

// " :::: x he i ,„ _ „ , 



det 



,i m 

/ d\\[(\ + g c )e-i- 

1 c=l 



(30) 



i.j — l.n 



with g c = 2/T c — 1 and 6ij — N(fa — fa) — i(yi + yj). The expression above coincides in every 
detail with that obtained in |q] for random GUE matrices deformed by a finite rank anti-Hcrmitian 



(27) 
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perturbation.^]. This completes the proof of universality for finite-rank deviations. 
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